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1 .  INTRODUCTION  - -  ~ ' 

recent  years,  increasing  interest  has  centered  on  the  development 
of  reliable  and  computationally  inexpensive  a  posteriori  error  estimates 
for  finite  element  computations,  (see  for  example  [1],  [2],  [3],  [4],  [5], 

[6] ,  [7])."^Such  estimates  can  provide  some,  often  critically  important, 
information  about  the  accuracy  and  reliability  of  the  computed  solution  as 
a  model  of  the  behavior  of  the  physical  phenomena  under  study.  At  the  same 
time  it  has  become  widely  accepted  that  these  estimates  also  constitute  a 
basic  tool  in  the  construction  of  efficient  adaptive  finite  element 
processes  which  are  designed  to  achieve  a  desired  error  tolerance  at  minimal 
cost  or  a  best  possible  solution  within  an  allowable  cost  range,  (see  e.g. 

[8],  [9L  [10],  [11],  [12]). 

/Up  to  now  most  of  this  work  concerned  linear  problems.  Not  unexpectedly, 
for  non-linear  problems  the  situation  is  much  more  difficult  and  the  theory 
is  by  far  not  as  well  developed.  In  part  this  is  due  to  the  many  special 
features  of  non-linear  problems  not  present  in  the  linear  case.  In  particular. 


such  problems  usually  involve  a  number  of  intrinsic  parameters  and  —  because 
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of  the  non-linear  nature  --  interest  centers  rarely  on  the  determination 
of  a  few  specific  solutions  for  fixed  parameter  values  but  instead  on  a 
more  general  study  of  these  solutions  under  various  changes  of  the  para¬ 
meters.  ^  In  structural  mechanics  the  parameters  may  characterize,  for 
instance,  load  points  and  load  directions,  material  properties,  geometrical 
data,  etc.,  and  the  set  of  all  solutions  depending  on  these  parameters  has 
been  called  the  equilibrium  surface  of  the  structure,  [13].  This  equili¬ 
brium  surface  provides  considerable  insight  into  the  behavior  of  the 
structure  and  its  stability  properties,  and  the  computational  task  is  to 
analyze  its  shape  and  characteristic  features.  For  this  the  principal 
tools  are  the  continuation  methods,  or  incremental  methods  as  they  are 
often  called  in  the  engineering  literature.  Generally,  these  processes 
allow  for  the  trace  of  any  path  on  the  surface  defined  by  a  parameter 
combination  with  one  degree  of  freedom. 

In  the  analysis  of  a  structural  problem  by  the  finite  element  method 
we  are  able  to  compute  only  approximate  points  along  a  path  on  the  solution 
surface  of  some  discretized  form  of  the  original  problem.  Then  we  are 
faced  with  the  need  for  assessing  and  controlling  the  errors  along  an  entire 
segment  of  such  a  path.  In  [14],  and  [15]  it  has  been  shown  that  effective, 
a  posteriori  error  estimates  and  adaptive  procedures  can  be  constructed 
which  meet  these  aims  and  which  can  be  successfully  incorporated  into  a 
general  continuation  process  for  tracing  paths  on  the  equilibrium  surface. 
But  these  results  were  essentially  restricted  to  problems  in  one  space 
dimension,  primarily  because  the  estimates  used  there  were  computationally 


relatively  expensive. 


Y*n  this  paper ^-we-  present^a  new  approach  to 


the  construction  of 
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a  posteriori  error  estimates  for  non-linear  problems  which  is  highly 
effective  and  at  the  same  time  computationally  rather  inexpensive. v  In 
fact,  these  estimates  use  a  linearized  form  of  the  problem  and  hence  their 
computation  can  be  accomplished  about  as  rapidly  as  in  the  case  of  linear 
problems.  Moreover,  the  approach  allows  us  to  bring  to  bear  most  of  the 
earlier  cited  results  about  linear  a  posteriori  error  estimates  and  hence 
it  applies  also  to  problems  in  more  than  one  space  dimension.  Based  on 
these  estimates,  a  prototype  software  system  for  the  adaptive  finite  element 
solution  of  a  class  of  two-dimensional,  parametrized  non-linear  problems  is 
now  under  construction  at  the  University  of  Pittsburgh.  It  was  dubbed 
NFEARS,  short  for  Nonlinear  Finite  Element  Adaptive  Research  Solver,  in 
analogy  to  an  earlier  developed  system  of  the  same  type  for  linear  problems 
called  FEARS,  (see  [8],  [9],  [11],  [16],  [17]). 

In  Section  2  below  we  present  the  general  principles  behind  the  new 
error  estimates.  Then  Section  3  outlines  some  of  the  features  of  the 
design  of  NFEARS.  Finally  in  Section  4  we  give  some  numerical  results  which 
show  the  effectiveness  of  the  error  estimates  and  of  the  adaptive  approach. 


2.  A  POSTERIORI  ERROR  ESTIMATES 

As  mentioned  in  the  Introduction,  the  new  error  estimates  for  non¬ 
linear  problems  utilize  the  earlier  developed  estimation  theory  for  linear 
problems.  It  would  be  impossible  to  present  here  an  account  of  that  theory; 
but  it  may  be  useful  to  illustrate  the  ideas  on  a  simple  example. 

Consider  the  two-point  boundary  value  problem 


L[u]  i  -  —  a($)  |^-+  b(s)u  =  c(s),  0  <  s  <  1, 


(2.1) 
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u(0)  =  u(l)  =  0  (2.2) 

where  the  coefficient  functions  a,  da/ds,  b,  and  c  are  assumed  to 
be  continuous  on  I  =  [0,1]  and  such  that 

0  <  aQ  ^  a(s)  <_  a-|  <  °°»  0  <_  b ( s )  <_  o ^  <  °°>  V  s  e  I.  (2.3) 

It  is  well  known  that  the  bilinear  form 

B(u,v)  =  f  [a(s)u‘v*  +  b(s)u,v]ds  (2.4) 

J 
I 

is  defined  and  continuous  on  H^(I)  *  H^(I),  where  in  (2.4)  primes  denote 
derivatives  with  respect  to  s.  Then  the  energy  norm  is  given  by 

||v||E  =  B(v,v)1/2.  (2.5) 

The  weak  solution  of  our  problem  is  now  the  unique  function 
uo  e  f°r  which 

B(uq,v)  *  F(v)  5  fvds,  VveHj(I).  (2.6) 

I 

Suppose  that  we  use  piecewise  linear  elements  on  some  mesh 


A:  0=sQ  <  s^  <  Sg  <  ...  sn+j  *  1,  n  =  n(A), 


(2.7) 


with  not  necessarily  uniform  steps  hfc  =  sfc  -  skl,  k  =  l,...,n+l.  In 
other  words  we  restrict  consideration  to  the  finite  dimensional  subspace 


S(A)  -  (u  e  hJ(I),  u(s) 


n 

l  x,  *|(s);  SE  I) 
1=1  1  1 


(2.8) 
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where  for  each  1  *  l,...,n,  <J>.  denotes  the  continuous  piecewise  linear 
"hat  function"  which  is  1  for  s  =  s^  and  zero  at  all  other  nodes  of 
A.  Then  the  finite  element  solution  uQ  e  S(a)  is  uniquely  defined  by 
the  condition 


B(uq,v)  =  F(v),  V  v  e  S(A). 

In  order  to  estimate  the  norm  of  the  finite  element  error  e  =  u  -  u 

o  o 

let  P.  denote  the  orthogonal  projection  of  H^(I)  onto  the  subspace 

K  0 

{v  G  Hq(I);  v ( s )  =0,  V  s  $  (sk_rsk)} 

with  respect  to  the  scalar  product  defined  by  B(u,v).  Then  with  the 
error  indicators 


we  can  define  the  error  estimator 


(2.9) 


n+1 

e  -  u 

k=l 


2 /2 
V 


for  which  it  can  be  shown  that  e  <  1 1 e 1 1 £  <  Ce  with  a  constant  C  <  00 
that  depends  only  on  the  bounds  aQ,  o^,  in  (2.3)  but  not  on  a,b,c, 
and  a.  The  restriction  of  the  function  Pke  to  Ik  =  [sk1 .sk]  is  the 
unique  solution  wk  of  the  original  problem  on  Ik  for  which 

wk(sk-l}  “  wk(sk)  *  °- 

Obviously,  we  are  Interested  in  simple  approximations  of  the  wk. 


-• . 
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For  instance,  we  may  use  the  quadratic  finite  element  approximation 


\M  =  Pkzk(s).  zk^s^  =  T  ^s-sk  l^sk“s^’  s  £  !k* 

hk 


Then  we  obtain  the  approximate  error  indicator 


[a(s)u^z£  +  b(s)uQzk  -  c(s)zk]ds}‘: 


nk  =  NWkl 


V 


[a(s)zkzk  +  b(s)zkzk]ds 
{  f  Vk  ds)2 


V 


^  |  [a(s)zkzk  +  b(s)zkzk]ds 


(2.10) 


where  rk(s)  =  L[uq](s)  -  c(s),  s  e  Ik>  is  the  residual  of  uQ  on  Ik> 

Evidently,  there  are  various  other  forms  of  the  error  indicators  and 
hence  of  the  a  posteriori  error  estimate.  For  more  details  and  proofs  we 
refer  to  [  1  ],  [  3  ],  [4  ]. 

In  order  to  summarize  the  situation  in  the  two-dimensional  case, 
consider  the  problem  defined  by 

[  [(Vv)TM(s,t)Vu  -  c(s,t)uv]ds  dt  =  0,  V  v  e  H][(n),  (2.11) 

0 

where  M(s,t)  is  a  symnetric,  positive  definitive,  2x2-matrix  for  all 
(s,t)  e  n  *  [0,1]  x  [0,1]  and  Vu  ■  (u$*ut)T  Is  the  gradient. 

Suppose  that  we  use  bilinear  elements  on  a  uniform  mesh  on  5  with 
step  h  >  0.  Then  the  evaluation  of  the  residuals  of  the  finite  element 
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solution  u„  results  in  a  Dirac  function  with  linear  distribution  on 
o 

the  sides  of  each  element.  More  specifically,  let  P  be  any  interior 
node  of  the  mesh  and  denote  its  four  horizontal  and  vertical  neighbors 
by  W,N,E,S  in  the  sense  of  the  windrose.  If  the  corresponding  function 
values  are  uQ(P),  uQ ( W ) ,  uQ (N) ,  uq ( E) ,  uQ(S),  respectively,  then  the 
following  jump  values  may  be  computed 


JH(P)  =  l  (0o(E)-Go(P))  -  1  (G0(P)-Q0(W)), 


JV(P)  =  l  (Go(N)-Go(P))  -  1  (u0(P)-u0(S)). 


At  any  point  P  on  one  of  the  vertical  boundaries  but  not  at  a  corner 
point  of  the  domain,  let  P'  be  the  immediate  horizontal  neighbor  of  P. 
Then  we  set  JH(P)  =  JH(P')  and  JV(P)  =  0.  Similarly,  for  any  point  P 
on  a  horizontal  boundary  but  not  at  a  corner  of  the  domain  we  set 
0V(P)  =  JV(P')  and  JH(P)  =  0  where  now  P'  is  the  immediate  vertical 
neighbor  of  P.  With  these  quantities  the  first  part  of  the  error  indicator 
of  any  element  t  of  our  mesh  is  given  by 


nl(x)  =  h  a  l  CM1l(s*t)JH(P)Z  +  M22(s,t)JV(P)Z]. 


(2.12) 


Here  M. ^  denote  the  elements  of  the  matrix  M,  the  sum  extends  over  all 
nodes  P  of  x  which  are  not  comers  of  n,  and  the  factor  a  equals 
4/3  if  t  has  a  node  at  a  corner  of  0  and  1,  otherwise. 

For  the  computation  of  the  residual  contribution  to  the  error  indicator 
let  R  be  given  by 
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R  =  *-3s  M12  +  at  M22^3t  u0  +  ^as  Mll  +  at  M2l  ^as  uo 


3  r. 


+  [m12  +  m91  ] 


21 Jasat  uo  '  c  uo 


where  the  M. .  are  again  the  elements  of  the  matrix  H.  Moreover,  let 

'  J 

R.  denote  the  values  of  R  at  the  four  Legendre-Gauss  points  of  the 

J 

element  t.  Then  the  residual  contribution  to  the  error  indicator  of  t 
is 


.4  4 


n2(x)  -  \  l  [(R,  -  i  (  l  R,)]* 

1  4/  j=l  J  4  j-1  J 


The  total  error  indicator  for  the  element  t  is  now 


n(-r)  =  n^x)  +  n2(x) 

and  the  sum  of  all  the  error  indicators  over  the  elements  of  the  mesh  is  the 
square  of  the  desired  error  estimate  e.  For  proofs  we  refer  to  [6]. 

In  [14],  and  [15]  the  above  indicated  approach  was  applied  directly 
to  nonlinear  problems.  In  that  case,  the  auxiliary  problems  for  the  w^ 
become  nonlinear  and  this  is  the  source  for  the  computational  expense 
mentioned  in  the  Introduction. 

In  order  to  reduce  this  expense  we  utilize  here  a  basic  property  of 
the  derivative  of  our  nonlinear  operators.  In  general,  the  parameterized 
equations  under  study  have  the  form 

F(y,X)  *  0 


(2.13) 


where  y  represents  a  state  variable  varying  in  an  appropriate  normed 
space,  and  X  is  a  finite  dimensional  parameter  variable.  For  simplicity 
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we  consider  here  only  the  case  of  a  non-singular  point  (u0>*0)  where 
the  derivative  DuF(uq »^0)  of  F  with  respect  to  the  state  variable  has 
a  continuous  inverse.  At  such  a  point  the  linearized  operator  has  the 
form 


L<VVM  =  F(Vxo}  +  V^VV^o5* 

Let  (uQ.Xo)  be  the  exact  solution  of  F(uo,AQ)  =  0  and  wQ  the 
solution  of  the  linear  equation  L(uQ ,XQ ) [wj  =  0.  Then  under  s  .a^le 
conditions  for  F  it  follows  that 

lluo-w0N  -  °(l|u0*i3iol,)  as  ^  *uo_Go  1 1  “"  °*  (2.14) 

In  our  setting,  uQ  represents  the  computed  finite  element  solution  and 
uQ  the  exact  solution  of  the  given  nonlinear  problem.  Now  (2.14)  implies 
that 


IIVa0N  =  +  o(1>>  as  HvDoll  (2.15) 

and  hence,  the  error  I |w  -ul I  between  the  solution  of  the  linearized- 
problem  and  the  computed  solution  is  asymptotically  equal  to  the  desired 
error  ||uo-uo||.  For  the  approximation  of  | I wQ-uo | |  we  may  apply  the 
earlier  developed  a  posteriori  error  estimates  for  linear  problems. 

In  order  to  illustrate  the  approach  we  consider  the  following  non¬ 
linear  version  of  the  problem  (2.1 )-(2.2) : 

N[u]  =  -  A(§|)  +  B(s,u)  -  C(s,A),  0  <  s  <  1 ,  (2.16) 


u(0)  =  u(l )  »  0. 


(2.17) 


A  weak  formulation  is 


[A(u')v'  +  B(s,u)v  -  C(s,A)v]ds  =0,  V  v  e  h](I) 
iJ  o 

where  O^B  denotes  the  derivative  of  B  with  respect  to  u.  Thus,  the 
linearized  problem  has  the  form 

[A'fu'JwV  +  DB(s,u)wv  -  C(s,X)v]ds  =  0,  V  v  e  hJ(I).  (2.18) 

t  U  0 

If  we  use  again  piecewise  linear  elements  on  the  mesh  (2.6),  then  we 
have  to  apply  the  linear  a  posteriori  estimates  to  the  linearized  problem 
(2.18).  More  specifically,  we  have  to  use  the  linear  problem  (2.4)  with 


a(s)  =  A'(u'(s)),  b(s)  =  DuB(s,u0),  c(s)  =  C(s,XQ). 


Hence,  if  we  proceed  as  in  (2.10)  then  we  obtain  the  error  indicators 


nk  = 


V 


f  [A1 (^q)UqZ^  +  B(s ,un)unzL  -  C(s,Xn) ]ds 


o'  o  k 


l"j  1  [A'(iio)zkzk  +  B(s,D0)zkzk]ds]|1 
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(2.19) 


where  the  denominators  are  assumed  to  be  non-zero.  A  more  complete  theory  of 
this  approach  will  be  given  elsewhere.  The  approach  is  exactly  the  same 
for  problems  involving  more  than  one  space  dimension. 


3.  THE  NFEARS  DESIGN 

As  mentioned  already  in  the  Introduction  a  prototype  software  system, 
called  NFEARS,  is  currently  under  development  at  the  University  of  Pittsburgh 


n 


which  has  the  following  characteristics: 

(i)  The  system  constitutes  an  applications- independent  finite 
element  solver  for  a  class  of  two-dimensional,  parametrized 
nonlinear  boundary  value  problems  defined  by  a  weak  variational 
formulation. 

( i i )  Adaptive  approaches  are  employed  throughout  and  the  a  posteriori 
error  estimates  outlined  in  Section  2  above  are  used  to  control 
the  process. 

(iii)  The  system's  design  is  analogous  to  that  of  the  linear  adaptive 
finite  element  solver  FEARS  described  in  [8],  [9],  [16]  and  [17]. 

Details  of  the  design  of  NFEARS  may  be  found  in  [18].  Accordingly  we  will 
outline  here  only  some  of  the  principal  features.  Moreover,  since  many 
design  aspects  of  NFEARS  correspond  to  those  of  FEARS  we  refer  also  to  the 
cited  references  for  that  system. 

The  permissible  domains  are  of  the  same  type  as  in  FEARS.  In  brief, 

2 

the  domain  Q  c.  R  is  the  union 

—  0*1  Li  &2  ^7  *  •  •  O 

2 

of  finitely  many  closed,  bounded  subsets  n.  c  R  with  disjoint,  non- 

J 

empty  interiors  Q-.  For  each  n.  a  one-to-one,  smooth  mapping  to.  onto 

J  J  J 

the  unit  square  [0,1]  *  [0,1  J*  of  R  is  given  and  these  mappings  satisfy 
certain  natural  compatibility  conditions  (see  e.g.  [9],  [16]  or  [17]). 

The  system  will  have  two  modules,  namely  for  the  solution  of  problems 
with  one  and  two  unknown  functions,  respectively.  In  the  case  of  one 
unknown  function  we  seek  a  function  u  defined  on  such  that 

(i)  u  satisfies  prescribed  boundary  conditions  on 

The  form  of  these  conditions  Is  analogous  to  those  of  FEARS 
and  Includes  inhomogeneous  Dirlchlet  conditions,  Neumann 
conditions,  as  well  as  mixed  boundary  conditions. 
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( i i )  u  is  a  stationary  point  of  the  functional 


I  A(I, ,I,,X,u)dsdt  -  2  f  f(s,t,A,u)u  dsdt 
ftJ  1  £  sii 


where 


(3.1) 


^1  =  ^  = 


are  invariants  with  respect  to  rotations  of  the  coordinate 
system. 


Similarly,  in  the  case  of  two  unknown  functions  we  are  seeking  a  vector 
u  =  (u-j.Up)  of  unknown  functions  on  ft  such  that  again  the  above 
conditions  (i)  and  (ii)  hold  but  with  the  functional  (3.1)  replaced  by 


f  A(I, , I9,I~,A,p)dsdt 
ftJ  12  3 


2 

ftJ 


C^i  (s,t,A,u)u<| 


(3.2) 


*  Cgfs.t.A.yJUgjdsdt 


where  now 


I'j  -  a^j  (vu-j  jVUg ) ,  Ig  =  agWu-j  .VUg) ,  1^  =  a^tu). 

This  class  of  problems  is  fairly  general  and  includes,  in  particular, 
most  of  the  basic  problems  of  elasticity  theory.  Two  parameters  A  and 
v  are  incorporated  in  the  formulation  and  hence  the  equilibrium  surface 
of  the  problems  under  study  Is  two-dimensional. 

The  fundamental  computational  process  Is  a  continuation  process  which 
allows  for  a  trace  of  any  path  on  the  surface  specified  by  a  given  para¬ 
meter  combination  with  one  degree  of  freedom.  More  specifically,  a  form 


of  the  PITCON  system  described  in  [19]  and  [20]  is  used  for  this  purpose. 

The  designs  of  the  data  structure  for  the  meshes  and  of  the  corres¬ 
ponding  access  algorithms  follow  essentially  those  of  FEARS  except  that 
in  NFEARS  a  biquadratic  element  will  be  used.  Some  details  on  this  and  on 
further  aspects  of  the  NFEARS  design  are  given  in  [18]. 


4.  SOME  NUMERICAL  EXAMPLES 

As  a  first  example,  consider  the  nonlinear  boundary  value  problem 
(2.16)-(2.17)  with 


A(t)  =  t/(l+t),  B(s,t)  =  0,  C(s,A)  =  A 


(4.1) 


in  which  case  the  exact  solution  is 


u(s)  =  -s  +  y  ln[(eX-l)s  +1],  0<  s  <  1 


For  growing  A  this  solution  increases  rapidly  within  a  small  interval 
near  s  =  0. 

A  continuation  process  has  been  used  to  compute  the  solution  path 
with  initial  point  at  the  origin  for  a  =  0.  The  points  with  parameter 
values  a  *  1.5  and  a  *  3.0  were  used  as  target  points.  At  these  points. 
Table  1  gives  a  comparison  of  the  exact  error  and  the  computed  a  posteriori 
error  estimates.  Uniform  meshes  with  degrees  of  freedom  n  *  4,8,16  were 
used  as  well  as  a  nonuniform  mesh  obtained  with  the  adaptive  procedure 
sketched  further  below.  The  table  shows  that  even  for  comparatively 
large  relative  errors,  the  effectlvlty  of  the  estimates  Is  excellent.  This 
has  been  our  general  experience  also  with  a  number  of  other  problems 
considered  so  far. 
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As  a  second  example  we  consider  the  two-dimensional  problem 

‘  ^A1(VU>  •  =  C(s,t,A),  V  (s,t)  e  a  =  (0,1)  x  (0,1)  (4.2) 

subject  to  the  boundary  conditions 

u  =  0  on  an.  (4.3) 

Evidently,  In  weak  form  the  linearized  problem  then  has  the  form  (2.11) 
with 


(D-jA-j  (Vu)  ^(VuA 

♦  C(s,t)  =  0. 

D/^Vu)  D2A2(7u )/ 

Here  D.jAj,  1,j  ■  1,2,  Is  the  derivative  of  the  coefficient  function  Aj 
with  respect  to  its  Ith-variable.  Thus  the  error  indicators  and  estimates 
can  be  computed  exactly  as  stated  In  Section  2. 

As  a  model  case  we  use  the  coefficient  functions 

A1  -  (us+ut)aV  A2  =  (us+ut)0V  {4-4) 

and  choose  C  in  such  a  way  that  the  exact  solution  of  (4. 2)- (4. 3)  is 

u  a  A  s(l-s)t(l-t),  V  s,t  e  ii.  (4.5) 

The  continuation  process  was  used  to  compute  the  solution  path  with  initial 
point  at  the  origin  for  A  *  0.  The  points  with  parameter  values  A  =  2,4.8,16 
were  used  as  target  points.  Uniform  meshes  with  m  *  4  and  m  =  16  elements 
were  used.  Table  2  presents  again  a  comparison  of  the  exact  error  and  the 


£ 


I  !e|  l/|  |u|  |%  e/MullX  e/|  |e|  |% 


A  =  2 

16 

64 

(|u||  =  0.2438 

0.5191  (-1) 

0.2551 (-1) 

0.5665(-l ) 

0. 261 5(-l ) 

21.29 

10.46 

23.23 

10.73 

109.1 

102.5 

A  =  4 

||u||  =  0.6896 

16 

0.1468 

0.1625 

21.29 

23.57 

110.7 

64 

0. 721 6(-l ) 

0.7449(-l ) 

10.46 

10.80 

103.2 

A  =  8 

|(u(|  =  1.950 

.!  16 

■i 

0.4153 

0.4726 

21.30 

24.23 

113.8 

64 

0.2041 

0.2137 

10.47 

10.96 

104.7 

A  =  16 

!|u||  =  5.517 

16 

1.175 

1.406 

21.29 

25.49 

119.7 

64 

0.5773 

0.6211 

10.46 

11.26 

107.6 

Table  2 

I 

computed  a  posteriori  error  estimates.  As  before,  the  table  reflects  our 
experience  that  the  effectivity  of  the  estimates  is  excellent  also  in  the 
two-dimensional  case  even  when  the  errors  are  large. 

It  is  now  widely  accepted  that  for  realistic  problems  it  is  rarely 
feasible  to  construct  numerical  processes  which  reliably  and  effectively 
achieve  a  desired  accuracy  at  reasonable  cost  and  yet  which  do  not  utilize 
some  form  of  adaptivity.  The  design  of  such  an  adaptive  procedure  depends 
strongly  on  the  goal  of  the  computation  (see  eg.  [7]);  but  in  all  cases 
the  availability  of  reliable  error  estimators  appears  to  be  central  to  the 
design  of  effective  adaptive  processes. 

In  nonlinear  problems  the  goal  of  the  computation  may  take  many 
different  forms.  For  example,  it  may  be  required  that  at  each  one  of  the 
computed  points  along  the  solution  path  the  error  estimate  does  not  exceed 
a  specified  tolerance.  In  other  cases,  the  goal  may  be  to  meet  the  prescriped 
error  tolerance  only  at  particular  target  points.  In  yet  other  cases, 
interest  may  center  on  the  accurate  calculation  of  certain  critical  points, 
such  as  buckling  points.  Alternately,  instead  of  the  computation  of  solution 
points  with  prescribed  error  behavior  we  may  focus  only  on  the  accuracy  of 
the  values  of  a  specified  functional  of  these  solutions,  such  as,  some 
stress  values,  etc. 

In  nonlinear  finite  element  computations,  the  principal  mechanisms  for 
the  control  of  the  adaptive  process  are  the  following: 

(1 )  Path  Controls 

(a)  Steplength  selection 

(b)  Local  parameter  selection 

(c)  Corrector  adjustments  (especially  In  the  case  of  iterative 
correctors) 
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(2)  Approximation  Controls 

(a)  Mesh-refinement 

(b)  Mesh-modification 

(c)  Changes  of  the  order  of  the  element 

The  path  controls  are  usually  part  of  the  continuation  process;  they  affect 
only  the  quality  of  the  numerical  solution  of  the  current  discretized 
problem.  The  errors  of  these  solutions  with  respect  to  the  exact  solution 
of  the  underlying  mathematical  model  can  only  be  influenced  by  means  of 
the  approximation  controls.  But  there  is  nevertheless  a  strong  interaction 
between  the  two  sets  of  controls  which  appears  not  to  have  been  addressed 
as  yet. 

There  have  been  only  relatively  few  efforts  of  incorporating  one  or 
several  of  the  approximation  controls  into  nonlinear  finite  element  solvers. 
Since  each  one  of  these  three  controls  has  certain  advantages  and  dis¬ 
advantages,  some  combination  of  them  may  well  be  desirable.  But,  even  in 
the  linear  case  the  construction  of  an  effective  combination  of  these  controls 
is  as  yet  not  fully  understood.  Accordingly,  the  current  design  of  NFEARS 
follows  the  model  of  FEARS  and  utilizes  only  mesh-refinements  for  the  control 
of  the  approximation  errors. 

Nevertheless,  the  design  of  an  effective  adaptive  mesh  refinement 
strategy  for  maintaining  a  prescribed  error  tolerance  at  all  points  of  the 
solution  path  still  remains  a  research  problem.  We  sketch  here  one  such 
design  which  appears  to  be  promising.  For  simplicity,  the  discussion  will 
be  restricted  to  the  one-dimensional  problem  (2.16)-(2.17). 

A  mesh  partition  function  Is  any  continuous,  strictly  monotone- 
increasing  function  <j>:  [0,1]  [0,1]  with  <j»(0)  *  0,  #(1 )  *  1 .  Then, 

for  any  n  >  1,  the  solutions  s^  of  the  equations 
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■  J-,  i  *  0,1,. ...n  (4.4) 

are  unique  and  define  a  mesh  (2.7)  on  our  domain  I  =  [0,1]  which  will 
be  denoted  by  a(<M). 

For  linear  problems  (2.1 )-(2.3)  it  was  shown  in  [3]  that 


with 


*(s)  =  ~  fS  (a(t)u  (t)2)1/3dt,  s  e  I 
Yo  Jo  0 


(4.5) 


(4.6) 


defines  an  asymptotically  optimal  partition  and  that  the  corresponding 
errors  for  the  meshes  An  =  A(<J>,n)  are 


|e| 


3/2 

To 

AI  n 


(1  + 


°<w>- 


as  h 


max 


0, 


(4.7) 


where  hmflx  denotes  the  maximal  mesh  step  h^  of  An-  The  asymptotic 
optimality  means  that  for  any  other  partition  function  \|»  and  all  n  for 
which  the  maximal  step  h„v  of  A(«i>,n)  is  sufficiently  small,  the  corres- 
ponding  error  is  not  less  thatn  (4.7). 

Let  ni»...»nn+^  denote  the  error  indicators  (2.9)  of  the  finite 
element  solution  for  a  given  mesh  A.  Then  it  turns  out  that 


"fS  '  M  1  J  ♦  °<W>-  «  hmax  *  »• 

!i 

This  can  be  used  to  obtain  from  the  computed  error  indicators  an 
approximation  of  the  optimal  partition  function  (4. 5)-(4. 6).  Evidently,  in 


the  nonlinear  case  we  may  proceed  analogously  by  using  again  the  linearized 
problem  as  we  did  in  the  computation  of  the  error  indicators. 
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As  mentioned  before,  our  aim  is  to  keep  the  error  estimate  at  each 
computed  point  below  a  given  error  tolerance.  As  tolerance-  it  is  useful 
to  define  the  quantity 


tol 


abs 


+  6 


rel 


ll“0N 


where  ||u  ||  is  the  norm  of  the  computed  solution  and  6abs  and  6rel 
are  absolute  and  relative  error  constants.  For  6gbs  =  0  this  means  that 
the  relative  error  has  to  be  less  than  6rel ,  while  for  6rgl  =  0  the 
absolute  error  is  not  permitted  to  exceed  6abs. 

If  at  any  point  along  the  solution  path  the  computed  solution  uQ 
has  an  error  estimate  below  0  tol,  where  0,  0  <  0  <  1,  is  a  given  factor, 

then  the  process  continues  with  the  same  mesh.  Otherwise  the  approximate 
optimal  mesh  partition  function  $  is  computed  from  the  error  indicators 

of  u  and  with  it  the  "ideal"  mesh  size  h  as  the  smallest  integer  not 

o 

below  Yq/2/(/T7  0  tol). 

This  allows  us  to  compute  an  "ideal"  mesh  A  *  A($,n)  and  to  construct 
from  the  current  mesh  A  a  refined  (or  de- refined)  mesh  A'  which  approxi¬ 
mates  A  in  some  way.  There  are  many  ways  to  accomplish  this  approximation 
and  a  detailed  discussion  of  such  techniques  shall  be  given  elsewhere.  Here 
we  present  only  some  typical  results  with  an  adaptive  mesh-refinement 
procedure  designed  along  this  line.  The  sample  problem  (2.16)-(2.17)  with 
the  coefficients  (4.1)  was  used  again.  For  the  run  given  In  Table  the 
tolerance  was  computed  with  6ab$  *  0.01,  5^  *  0.1,  and  the  tolerance 

factor  0  *  0.75. 
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X 

Nused 

N ideal 

error 

estlm 

exact 

error 

tolerance 

decision 

0 

4 

- 

0 

0 

- 

0.1231 

4 

4 

0.891 3(-2) 

0.8913(-2) 

0.1 344(-l ) 

proceed 

0.3395 

4 

7 

0.2516H) 

0. 251 7 ( - 1 ) 

0.1949(-1) 

refine 

7 

7 

0.1 504 ( - 1 ) 

0.1 504 ( - 1 ) 

0.1971 (-1 ) 

proceed 

0.5551 

7 

9 

0.2440(-l) 

0.2441 (-1) 

0.2594(-l ) 

ref i ne 

11 

9 

0.1 553(-l ) 

0.1 554 (-1  ) 

0.2604(-l) 

proceed 

0.7655 

11 

10 

0.21 33 ( -1 ) 

0.21 34 ( -1 ) 

0.3323(-l ) 

proceed 

1.15841 

11 

11 

0.3393(-l) 

0.3394(-l) 

0.4406(-l) 

proceed 

2.1161 

11 

15 

0.9654(-l) 

0.9701 (-1) 

0.7504(-l) 

refine 

15 

13 

0.5532(-l) 

0.5532(-l ) 

0. 7570( -1  ) 

proceed 

3.1192 

15 

21 

0.1603 

0.1620 

0.1146 

refine 

18 

15 

0.91 75 (- 1 ) 

0.91 67 (-1 ) 

0.1159 

proceed 

Table  3 


i 

* 

i 

a 

4 


Clearly,  as  expected,  whenever  a  point  Is  accepted  and  the  process  continues 
the  error  estimate  Is  below  the  prescribed  tolerance. 
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